Molecular dynamics of flows in the Knudsen regime 
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Novel technological applications often involve fluid flows in the Knudsen regime in 
which the mean free path is comparable to the system size. We use molecular dy- 
namics simulations to study the transition between the dilute gas and the dense 
fluid regimes as the fluid density is increased. 
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In a landmark paper |I] , Maxwell applied his kinetic theory |2j to the study of the slip length 
associated with fluid flow close to a solid surface. He postulated a linear combination of two 
kinds of behaviour when a molecule was incident on a wall - it is either specularly reflected or 
it undergoes enough collisions with the wall molecules so that it loses memory of the incident 
velocity and emerges with a velocity characteristic of the wall temperature. We revisit this 
problem with molecular dynamics simulations - as the properties of the wall are varied, we 
find both kinds of behaviour suggested by Maxwell, but the crossover from one to the other 
is not as simple as envisioned by him. In addition, rich behaviour is found in the fluid flow 
properties on varying the fluid density from the dense liquid regime to the Knudsen regime 
(large mean free path compared to the system size). Our studies of flows in the subcontinuum 
regime provide a unique window on the nature of boundary conditions at a fluid - solid 
interface and are relevant to the behaviour of nanoscale mechanical devices, high altitude 
aerodynamics, aerosol dynamics, isotope separation schemes, flow cooling techniques, and 
energy transfer in molecular collisions [0-0. 

Studies of flows at nanoscales involve questions pertaining to a sub-continuum regime [|7||| 
for which experiments are often unable to provide the requisite answers. The behaviour of 
materials at the nanoscale is often quite different from that in the bulk |9|-12|. Numerical 



simulations based on molecular dynamics (MD) offer a natural tool to handle the interacting 
many-body nanoscale physics. A detailed understanding of the behaviour of materials at very 
small scales is necessary for the fabrication of miniature devices and for the manipulation 
of materials at the molecular scale. Such microscopic studies are useful for assessing the 
range of validity of traditional continuum theories. In addition, atomistic simulations are an 
invaluable tool for deducing constitutive relationships, which may then be fed into continuum 
calculations. 

Here, we report on MD studies of the boundary conditions in single-fluid flows in the Knudsen 
regime in which the mean free path, A, becomes comparable to the characteristic device 
dimension, L. Such flows behave very differently from bulk, continuum expectations. There 
is considerable work on the modified boundary conditions based on either semi-heuristic 



kinetic theory arguments or using an analysis of the Boltzmann equation, but inevitably 



severe simplifications are made concerning the dynamics of the system ||14|| . Furthermore 



in such treatments, the solid boundaries are not handled in a realistic manner as having 



atomic structure on length scales comparable to that of the fluid constitutents [|T^,|nj • 
We study Lennard- Jones fluids undergoing two kinds of flow, Couette and Poiseuille, between 
two solid walls with an atomic structure. Three kinds of walls, two attractive and one 
repulsive, are discussed. We consider a range of fluid densities to study the crossover from 
flows described by continuum fluid mechanics to dilute gas flows in the Knudsen regime. 
In a Lennard- Jones fluid, two atoms separated by a distance r interact with the potential 



V LJ (r) = 4e 
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where a is the size of the repulsive core [the potential is truncated at 2.2a}. The flows 
take place between two solid planar walls parallel to the xy plane with periodic boundary 
conditions imposed along the x and y directions. Couette flow was generated by moving 
the walls at constant velocity in opposite directions along the x-axis whereas Poiseuille flow 
was induced by an application of uniform acceleration g on all the fluid molecules in the 
x-direction. 



Following Thompson and Robbins | I7| , the wall of atoms form two [001] planes of an fee 
lattice with each wall atom tethered to a fixed lattice site by a harmonic spring of spring 
constant k. Our studies where done for fc=400 so that the mean-square displacement about 
the lattice sites did not exceed the Lindemann criterion for melting but was close to the value 
set by this criterion. The wall fluid interactions were modelled by another Lennard- Jones 
potential, 



V wf (r) = 16e 
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with the parameter A varying between 1 and 0, corresponding to attractive and repulsive 
walls respectively. For the narrowest channel that we studied, the fluid atoms were confined 
to a space measuring 13.6cr by 5.1a in the xy plane and 12.75a between the walls. The 



"basic" walls in our studies contain 96 wall atoms each (in each wall plane, the atomic 
periodicity is 0.85a). We have also considered the case of walls which are about 2.5 times 
denser and consist of 176 wall atoms each - we shall refer to these as "dense" walls. 
The equations of motion were integrated using a fifth-order Gear predictor-corrector algo- 
rithm jl8| with a step size of 0.005r, where r = (ma 2 /e) 1//2 and m is the atomic (wall and 
fluid) mass. The fluid and wall temperatures were fixed at ksT/e— 1.1. We present results 
obtained with the use of Langevin thermostat ]19| but the results obtained with the use of 
the Nose-Hoover thermostat |]20|pT| were virtually identical. The noise was applied in all di- 
rections during equilibration and only in the y-direction during the data-taking phase. The 
averages were obtained over runs which lasted at least 5000 r, and in the extreme dilution 
case, for up to 400000 r. An initial time of at least 400r was spent for equilibration. The 
spatial averaging was carried out in slabs of width a/4 along the z-axis. 
In the dilute gas limit, the Knudsen number, Kn , is defined as X/L, i.e., as pa" 2 'L)^ 1 , 
where A is the mean free path, p is the fluid number density and L is the characteristic width 
of the channel. A study of the density profile shows qualitatively different behaviours of the 
density profiles for the attractive and repulsive walls. In the former case, two adsorbed layers 
form near both walls whereas there is merely a depletion zone, whose thickness is of order 
a, alongside the repulsive walls. The first layer at the attractive walls is fairly periodic, as 
measured by the static structure factor, and the resulting local order corresponds to either 
a square or hexagonal lattice depending on whether the wall is basic or dense. L is defined 
as an effective channel width in which the flow takes place (in the attractive case this is the 
distance between centers of the second layers). In the center of the channel, the fluid density 
is essentially constant and this value is chosen as p in the definition of Kn. 
Figure 1 shows the Poiseuille flow velocity profiles for three different values of p (or Kn) 
for the attractive wall case. We have confirmed with studies of channels of three different 
widths and the two wall densities that the crossover value of p above which continuum fluid 
mechanics holds, px, is somewhat less than 0.19, corresponding to a Kn of order 0.1 for the 
channel of the smallest width. In this regime, the profiles are clearly parabolic and essentially 
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with no slip at the walls. In contrast, the velocity profile at the highest Kn that we were 
able to study (corresponding to just 4 particles in the interior of the channel) is almost that 
of plug flow with a huge slip length. For repulsive walls, the plug flow regime extends to 
higher densities because such walls provide little coupling to the fluid and correspondingly, 
Pk becomes of order 0.3. Furthermore, in the continuum regime, the velocity field exhibits 
substantial slip. A consistent behaviour is found for Couette flow with higher slip for the 
repulsive walls (compared to the attractive wall case) for p > p K and essentially zero fluid 
velocity at the walls at the highest Kn. 

We now proceed to study the effect of interactions between the fluid particles undergoing 
Poiseuille flow on increasing p within the channel. In Poiseuille flow, the velocity field is 
highest in the middle of the channel and is denoted by v max . At the lowest densities, one 
expects that v max ought to scale as L. This is because the driving force is effective in 
accelerating the fluid particles only during the time spent traversing the channel between 
collisions with the walls (and the associated adsorbed layers). This time simply scales as the 
channel width. On the other hand, at large p, corresponding to usual fluid behaviour, v max 
ought to scale as L 2 as given by the Navier-Stokes equation with no-slip boundary conditions. 
Figure 2 shows a plot of how v max interpolates between the low and high Knudsen number 
behaviours as the fluid density is varied. Strikingly, t> max shows a maximum around px with 
the peak being more pronounced for larger L's. We have confirmed that this qualitative 
behaviour is essentially independent of the nature of the walls. Physically, px corresponds 
to the situation when the fluid intermolecular interactions begin to dominate leading to a 
diminished influence of the wall in the interior of the fluid. 

Figure 3 shows a vivid illustration of the role played by the intermolecular interactions in 
leading to an enhancement of t> max even as the fluid density is moderately increased from 
the high Kn limit. For a repulsive wall (not shown), the particle is reflected approximately 
specularly and moves out of the vicinity of the wall immediately. The bottom panel shows 
qualitatively different behaviour in the (basic) attractive wall case. Here, again as suggested 
by Maxwell, a particle either thermalizes in the vicinity of the wall for a relatively long time 
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and then escapes from it or occasionally still undergoes specular reflection. On increasing the 
fluid density, one gets the situation illustrated in the middle panel in which intermolecular 
collisions keep the particles longer within the channel thus reducing the time of localization 
near the walls and effectively increasing the velocity field in the middle of the channel (Figure 
2). At still higher densities (~ Pk\ the intermolecular interactions become sufficiently large 
so that liquid-like behaviour sets in and t> max begins to decrease. 



A key quantity that characterizes the fluid solid interface is the collision kernel ||i4| , |22|1 - the 
probability density that a molecule striking the surface with a given velocity reemerges with 
a specific velocity after a given residence time. For the repulsive wall, we find that Maxwell's 
specular reflection scenario holds exceedingly well and the residence time is negligible. On 
the other hand, for the attractive walls the distribution of residence times develops a long 
tail. Strikingly, the velocity distribution of a particle after the collision with the wall is 
substantially independent of the incoming velocity. As an illustration, Figure 4 shows the 
case of particles with unit speed incident normal to the wall. The inset shows the probability 
distribution of residence times whereas the main figure shows the distributions of the normal 
(z) and tangential (x) velocities of the scattered molecules. The results are in excellent accord 
with Maxwell's thermal wall scenario in which the probability densities corresponding to v x 
and v 2 are fli yi3l,|23 



and 
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<t>z(v z ) = t-= v z exp( :— =-) (4) 

respectively. As expected for a thermal wall, the residence time is found to be uncorrelated 
with the outgoing velocity. 

How does the behaviour crossover from that of specular reflection to that of a thermal wall? 
Following Maxwell 1 , a common hypothesis is a linear combination of both behaviours with 
weights (1 — /) and / respectively. Such an intermediate case is realized, e.g., on consid- 
ering the scattering from a wall with the wall-fluid potential given by Eq. 2 with A— 1/4. 



The velocity distributions of the scattered particles clearly deviate from pure thermal be- 
haviour and the data might be interpreted as containing a specular component. However, 



the resulting distributions are more complex than the simple linear combination form |24 



The behaviour of the collision kernel is further enriched on considering the role of the in- 
teractions between the fluid molecules as the Knudsen number decreases to a value of 0.27 
corresponding to a density of p = 0.066 (figure 5). The influence of other fluid molecules 
leads to the distributions more closely approaching the thermal ones. The nature of the 
boundary conditions at a fluid-solid interface is thus complex but akin to those envisioned 
by Maxwell. 

We are indebted to Mark Robbins and Riina Tehver for many valuable discussions. This work 
was supported by KBN (Grant No. 2P03B-025-13), NASA, and the Petroleum Research 
Fund administered by the American Chemical Society. 



7 



REFERENCES 

[1] J. C. Maxwell, Phil. Trans. R. Soc. London. Ser. A 170, 231 (1867) 
[2] J. C. Maxwell, Phil. Mag. 19, 20 (1860) 

[3] E. P. Muntz, D. Weaver, and D. Campbell, eds., Rarified Gas Dynamics, (American 
Institute of Aeronautics and Astronautics, Washington, 1989) 

[4] C.-J. Kim, ed., Microelectromechanical systems (MEMS), (American Society of Mechan- 
ical Engineers, New York, 1997) 

[5] W. H. Marlow, ed., Aerosols: Aerosol Microphysics I - Particle Interactions, (Springer, 
New York, 1980) 

[6] W. Ehrfeld, Elements of Flow and Diffusion Processes in Separation Nozzles, (Springer, 
New York, 1983) 

[7] J. Koplik and J. R. Banavar, Computers in Physics 12, 424 (1998) 

[8] J. Koplik and J. R. Banavar, Annu. Rev. Fluid Mech. 27, 257-92 (1995) 

[9] S. Sugano, Microcluster Physics, (Springer, Berlin, 1991) 

[10] B. Bhushan, J. N. Israelachvili, and U. Landman, Nature 374, 607 (1995) 

[11] I. Bitsanis, T. K. Vanderlick, M. Tirrel, and H. T. Davis, J. Chem. Phys. 89, 3252 
(1988) 

[12] W. Loose and S. Hess, Rheol. Acta 28, 91 (1989) 

[13] E. H. Kennard, Kinetic Theory of Gases, (McGraw-Hill, New York, 1938) 

[14] C. Cercignani, The Boltzmann Equation and its Applications (Springer, Berlin, 1988) 

[15] D. K. Bhattacharya and G. C. Lie, Phys. Rev. A 43, 761 (1991) 

[16] D. L. Morris, L. Hannon, and A. L. Garcia, Phys. Rev. A 46, 5279 (1992) 

s 



[17] P. A. Thompson and M. O. Robbins, Phys. Rev. A 41, 6830 (1990) 

[18] See, e.g. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon, 
Oxford, 1987) 

[19] G. S. Grest and K. Kremer, Phys. Rev. A 33, 3628 (1986) 

[20] S. Nose, Mol. Phys. 52, 255 (1984) 

[21] W. G. Hoover, Phys. Rev. A 31, 1695 (1985) 

[22] Y. Matsumoto, N. Yamanishi and H. Shobatake, Proceedings of the 19th International 
Symposium on Rarefied Gas Dynamics, 995 (1994) 

[23] R. Tehver, F. Toigo, J. Koplik, and J. R. Banavar, Phys. Rev. E 57, R17 (1998) 

[24] The velocity distributions in this case are shown (figure 10) in a conference paper, M. 
Cieplak, J. Koplik, and J. R. Banavar, Physica A 274, 281 (1999). The complemen- 
tary case of a repulsive wall is considered in detail in this work and the results are 
benchmarked against our basic findings presented here. 



Figure captions 



FIG. 1. Velocity profiles of Poiseuille flow with g = 0.01cr/r 2 for the smallest width channel 
with basic attractive walls. The values of the interior fluid density and the corresponding 
Knudsen number are indicated. The filled circles denote the highest Kn case. The inner wall 
layers are at the edges of the figure. Effectively, one observes no slip boundary conditions 
for the two higher density cases (the slip length is around 1.7 'a measured with respect to the 
wall location, but note the presence of the adsorbed layers at the wall). 

FIG. 2. Plot of the fluid velocity in the middle of the channel in Poiseuille flow with 
g = O.Olr/cr 2 as a function of fluid density (and the Knudsen number in the inset). The plot 
is for three different channel widths. The Lq in the figure denotes the narrowest channel. For 
this case, both attractive walls were considered - the basic (hexagons) and higher (asterisks) 
wall densities. 2L and 4L denote channels that are effectively twice and four times as wide 
as the smallest channel. 

FIG. 3. Plot of the time dependence of the z-coordinate of a typical fluid particle un- 
dergoing Poiseuille flow in the narrowest channel. The horizontal dotted lines indicate the 
^-coordinates of the wall molecules. The solid horizontal line refers to the centers of the 
adsorbed first layers. 

FIG. 4. The probability density distributions for the z (normal) and x components of 
the outgoing velocity of a fluid molecule after collision with the wall and the associated 
adsorbed layer. The solid histograms are obtained based on 3000 scattering events. The 
dotted lines are the analytic predictions of Maxwell, eqs. 3 and 4, discretized with a bin size 
of 0.1 ct/tq as in the simulations. The inset shows the distribution of residence times. The 
residence time is measured as the time interval between a fluid molecule leaving a force-free 
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region (near the wall) and its reemergence. The peak of the curve (at a non-zero value of 
r) corresponds to a virtually instantaneous bounce back of the fluid molecule and accounts 
for the time needed to approach the wall and return. 
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